Load the data

Read data from well 3 - Spleen.

seu <- readRDS("data/exp11W3_init.rds")

seu[["percent.mt"]] <- PercentageFeatureSet(seu, pattern = "^mt-")
seu[["percent.rt"]] <- PercentageFeatureSet(seu, pattern = "^Rp[ls]")

# Remove TRAV, TRBV, mito and ribo genes: spleen
seu_annots <- seu@meta.data
seu_counts <- seu[!grepl('Tr[ab]v|mt-|Rp[ls]|MGP-', rownames(seu)),]

Filtering out nrow(seu)-nrow(seu_counts) genes:

data.table(rownames(seu)[which(grepl('Tr[ab]v|mt-|Rp[ls]|MGP-', rownames(seu)))]) %>% 
  kbl() %>% 
  kable_styling(font_size = 11,
                bootstrap_options = c("striped", "hover", 
                                      "condensed", "responsive")) %>% 
  scroll_box(width = "100%", height = "400px")

Clustering and automated annotation

seu <- CreateSeuratObject(counts = seu_counts[["RNA"]]@counts, 
                          project = "seu", 
                          min.cells = 3, 
                          min.features = 200, 
                          meta.data = seu_annots)

  seu <- NormalizeData(object = seu)
  seu <- ScaleData(seu, verbose = FALSE)
  seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = 1000, 
                                       verbose = FALSE)
  seu <- RunPCA(seu, npcs = 20, verbose = FALSE)
  seu <- RunUMAP(seu, reduction = "pca", dims = 1:20)
  seu <- FindNeighbors(seu, dims = 1:20)
  seu <- FindClusters(seu, resolution = 0.9)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 8827
## Number of edges: 312856
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8394
## Number of communities: 18
## Elapsed time: 0 seconds

counts <- GetAssayData(seu[["RNA"]])

seu_singler <- CreateSinglerObject(counts=counts,
  project.name= "E07_Tregs", # choose
  min.genes = 200, # ignore cells with fewer than 200 transcripts
  technology = "10x", # choose
  species = "Mouse",
  citation = "Our 10x data exp 01-05", # choose
  normalize.gene.length = FALSE,        # needed for full-length platforms (e.g. smartseq)
  variable.genes = "de",  # see vignette
  fine.tune = FALSE, # TRUE would take very long
  reduce.file.size = TRUE, # leave out less-often used fields 
  do.signatures = FALSE,
  do.main.types = TRUE,
  numCores = 16)

## [1] "Dimensions of counts data: 12934x8827"
## [1] "Annotating data with Immgen..."
## [1] "Variable genes method: de"
## [1] "Number of DE genes:2610"
## [1] "Number of cells: 8827"
## [1] "Number of DE genes:2610"
## [1] "Number of clusters: 10"
## [1] "Annotating data with Immgen (Main types)..."
## [1] "Number of DE genes:1773"
## [1] "Number of cells: 8827"
## [1] "Number of DE genes:1773"
## [1] "Number of clusters: 10"
## [1] "Annotating data with Mouse-RNAseq..."
## [1] "Variable genes method: de"
## [1] "Number of DE genes:2702"
## [1] "Number of cells: 8827"
## [1] "Number of DE genes:2702"
## [1] "Number of clusters: 10"
## [1] "Annotating data with Mouse-RNAseq (Main types)..."
## [1] "Number of DE genes:2134"
## [1] "Number of cells: 8827"
## [1] "Number of DE genes:2134"
## [1] "Number of clusters: 10"

seu <- AddMetaData(seu,
seu_singler$singler[[1]]$SingleR.single$labels,
                                 col.name = "Immgen_annot_single")
seu <- AddMetaData(seu,
seu_singler$singler[[1]]$SingleR.single.main$labels,
                                 col.name = "Immgen_annot_single_main")
seu <- AddMetaData(seu,
seu_singler$singler[[2]]$SingleR.single.main$labels,
                                 col.name = "MouseRNAseq_single_main")
seu <- AddMetaData(seu,
seu_singler$singler[[2]]$SingleR.single$labels,
                                 col.name = "MouseRNAseq_single")

cowplot::plot_grid(DimPlot(seu, group.by = "MouseRNAseq_single_main"),
    DimPlot(seu, group.by = "Immgen_annot_single_main"), ncol = 2)

seu_spleen <- seu

DimPlot(seu_spleen, group.by = "Immgen_annot_single_main", raster = T, raster.dpi = c(900,900))

#ggsave("./fig/spleen_annotations.svg", width = 16, height = 12, units = "cm")

seu_spleen_filt <- subset(seu_spleen, nFeature_RNA > 750 & percent.mt < 7.5 & seurat_clusters %in% 0:14 &
                          Immgen_annot_single_main %in% c("T cells", "NK cells", "NKT", "Tgd", "ILC")  )

seu_spleen$keep <- colnames(seu_spleen) %in% colnames(seu_spleen_filt)


DimPlot(seu_spleen, group.by = "keep", raster = T, raster.dpi = c(900,900), pt.size = 3.5, cols = c("red","grey88"))

#ggsave("./fig/spleen_qc.svg", width = 13, height = 12, units = "cm")

Add metadata

We need to connect the hashtags of the analyzed mice with MouseIDs and matching diets.

seu_spleen$hashtags %>% table

## .
## Mouse_Hashtag_1 Mouse_Hashtag_2 Mouse_Hashtag_3 Mouse_Hashtag_4 Mouse_Hashtag_5 
##             942             995            1004             936            1044 
## Mouse_Hashtag_6 Mouse_Hashtag_7 Mouse_Hashtag_8 Mouse_Hashtag_9 
##             990             988             949             979

md_spleen <- seu_spleen@meta.data 
md_spleen <- md_spleen %>% mutate(MouseID = str_replace_all(string = hashtags, pattern = "_Hashtag", replacement = ""))
md_spleen <- md_spleen %>% mutate(Diet = if_else(MouseID %in% c("Mouse_1", "Mouse_2", "Mouse_3", "Mouse_4", "Mouse_5"), "GFD", "Standard"))

seu_spleen@meta.data <- md_spleen

DimPlot(seu_spleen, label = T)

Remove contaminating and low quality cells

VlnPlot(seu_spleen, features = c("nCount_RNA","nFeature_RNA","percent.mt","percent.rt"), ncol = 2)

seu_spleen_filt <- subset(seu_spleen, nFeature_RNA > 750 & percent.mt < 7.5 & seurat_clusters %in% 0:14)

Filtering removed ncol(seu_spleen)-ncol(seu_spleen_filt) cells

seu_spleen_filt <- NormalizeData(object = seu_spleen_filt)
  seu_spleen_filt <- ScaleData(seu_spleen_filt, verbose = FALSE)
  seu_spleen_filt <- FindVariableFeatures(seu_spleen_filt, selection.method = "vst", nfeatures = 500, 
                                       verbose = FALSE)
  seu_spleen_filt <- RunPCA(seu_spleen_filt, npcs = 12, verbose = FALSE)
  seu_spleen_filt <- RunUMAP(seu_spleen_filt, reduction = "pca", dims = 1:12)
  seu_spleen_filt <- FindNeighbors(seu_spleen_filt, dims = 1:12)
  seu_spleen_filt <- FindClusters(seu_spleen_filt, resolution = 0.7)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 8236
## Number of edges: 269011
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8524
## Number of communities: 13
## Elapsed time: 0 seconds

  DimPlot(seu_spleen_filt, group.by = "Immgen_annot_single_main", label = T, repel = T) 

  DimPlot(seu_spleen_filt, group.by = "MouseRNAseq_single_main", label = T, repel = T) 

  DimPlot(seu_spleen_filt, label = T, repel = T) 

 seu_spleen_filt <- reduced_annot(seu_spleen_filt, "Immgen_annot_single", n_cells_to_other = 30)

DimPlot(seu_spleen_filt, group.by = "Immgen_annot_single_red", label = T, repel = T) 

  ##ggsave("fig/Immgen_dimplot_spleen_filt_Immgen_single.png", width = 40, height = 20, units = "cm")

Analysis of spleen clusters

Immgen Annotations

seurat_meta_data <- seu_spleen_filt@meta.data

# Create grouped dataframe, calculate the frequencies of clusters
df4 <- seurat_meta_data %>% 
  group_by(MouseID, Immgen_annot_single_red) %>% 
  summarise(n = n()) %>% 
  unique() %>% 
  mutate(freq = n / sum(n)) 

# Control - all sums should be one
# df4 %>% group_by(MouseID) %>% summarise(sum = sum(freq))

# As we've lost non-grouping variables, let's join them back
md_to_join <- seurat_meta_data %>% 
  dplyr::select(MouseID, Diet) %>% 
  unique()

df4  <- left_join(df4, md_to_join)

# The final plot

df4 %>% 
  ggplot(aes(x = Diet, y = freq)) + # you can change the x to whatever variable you're interested in
  geom_boxplot(outlier.shape = NA) + 
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0) + 
  geom_jitter(position=position_jitter(0.2), size = 2, aes(color = Diet)) + # in aes, you can also use shape or fill (for the shapes that allow it)
  facet_wrap(~factor(Immgen_annot_single_red, labels = str_replace(levels(factor(df4$Immgen_annot_single_red)), pattern = " \\(", replacement = " \n(")), scales = "free", ncol = 5) +
  ylab("Frequency") +
  xlab("Condition") +
  theme_classic() +
  theme(strip.background = element_blank(), panel.grid = element_blank()) + 
  theme(axis.text.x = element_text(angle = 90)) + ggtheme() 

 ##ggsave("fig/clusters_boxplots_spleen_filt_Immgen_single.png", width = 27, height = 35, units = "cm")

Seurat clusters

seurat_meta_data <- seu_spleen_filt@meta.data

# Create grouped dataframe, calculate the frequencies of clusters
df4 <- seurat_meta_data %>% 
  group_by(MouseID, seurat_clusters) %>% 
  summarise(n = n()) %>% 
  unique() %>% 
  mutate(freq = n / sum(n)) 

# Control - all sums should be one
# df4 %>% group_by(MouseID) %>% summarise(sum = sum(freq))

# As we've lost non-grouping variables, let's join them back
md_to_join <- seurat_meta_data %>% 
  dplyr::select(MouseID, Diet) %>% 
  unique()

df4  <- left_join(df4, md_to_join)

# The final plot

df4 %>% 
  ggplot(aes(x = Diet, y = freq)) + # you can change the x to whatever variable you're interested in
  geom_boxplot(outlier.shape = NA) + 
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0) + 
  geom_jitter(position=position_jitter(0.2), size = 2, aes(color = Diet)) + 
  facet_wrap(~seurat_clusters, scales = "free", ncol = 5) +
  ylab("Frequency") +
  xlab("Condition") +
  theme_classic() +
  theme(strip.background = element_blank(), panel.grid = element_blank()) + 
  theme(axis.text.x = element_text(angle = 90))

 ##ggsave("fig/clusters_boxplots_spleen_filt_clusters.png", width = 16, height = 16, units = "cm")

Filter lymphoid populations

Filter high quality & lymphoid populations

seu_spleen_filt_lymph <- subset(seu_spleen, nFeature_RNA > 750 & percent.mt < 7.5 & Immgen_annot_single_main %in% c("T cells", "NK cells", "NKT", "Tgd", "ILC"))

seu_spleen_filt_lymph <- NormalizeData(object = seu_spleen_filt_lymph)
  seu_spleen_filt_lymph <- ScaleData(seu_spleen_filt_lymph, verbose = FALSE)
  seu_spleen_filt_lymph <- FindVariableFeatures(seu_spleen_filt_lymph, selection.method = "vst", nfeatures = 500, 
                                       verbose = FALSE)
  seu_spleen_filt_lymph <- RunPCA(seu_spleen_filt_lymph, npcs = 12, verbose = FALSE)
  seu_spleen_filt_lymph <- RunUMAP(seu_spleen_filt_lymph, reduction = "pca", dims = 1:12)
  seu_spleen_filt_lymph <- FindNeighbors(seu_spleen_filt_lymph, dims = 1:12)
  seu_spleen_filt_lymph <- FindClusters(seu_spleen_filt_lymph, resolution = 0.8)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 7517
## Number of edges: 245163
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8281
## Number of communities: 15
## Elapsed time: 0 seconds

  DimPlot(seu_spleen_filt_lymph, group.by = "Immgen_annot_single_main", label = T, repel = T) 

  DimPlot(seu_spleen_filt_lymph, group.by = "MouseRNAseq_single_main", label = T, repel = T) 

  DimPlot(seu_spleen_filt_lymph, label = T, repel = T) 

Analysis of spleen clusters

Annotations by diet

seurat_meta_data <- seu_spleen_filt_lymph@meta.data

# Create grouped dataframe, calculate the frequencies of clusters
df4 <- seurat_meta_data %>% 
  group_by(MouseID, seurat_clusters) %>% 
  summarise(n = n()) %>% 
  unique() %>% 
  mutate(freq = n / sum(n)) 

# Control - all sums should be one
# df4 %>% group_by(MouseID) %>% summarise(sum = sum(freq))

# As we've lost non-grouping variables, let's join them back
md_to_join <- seurat_meta_data %>% 
  dplyr::select(MouseID, Diet) %>% 
  unique()

df4  <- left_join(df4, md_to_join)

# The final plot

df4 %>% 
  ggplot(aes(x = Diet, y = freq)) + # you can change the x to whatever variable you're interested in
  geom_boxplot(outlier.shape = NA) + 
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0) + 
  geom_jitter(position=position_jitter(0.2), size = 2, aes(color = Diet)) + # in aes, you can also use shape or fill (for the shapes that allow it)
  facet_wrap(~seurat_clusters, scales = "free", ncol = 5) +
  ylab("Frequency") +
  xlab("Condition") +
  theme_classic() +
  theme(strip.background = element_blank(), panel.grid = element_blank()) + 
  theme(axis.text.x = element_text(angle = 90))

##ggsave("fig/clusters_boxplots_spleen_filt_lymph_clusters.png", width = 20, height = 15, units = "cm")

Markers spleen

markers_spleen <- FindAllMarkers(seu_spleen_filt_lymph, logfc.threshold = log(2))

markers_spleen %>% 
  kbl() %>% 
  kable_styling(font_size = 11,
                bootstrap_options = c("striped", "hover", 
                                      "condensed", "responsive")) %>% 
  scroll_box(width = "100%", height = "400px")
#write.csv(markers_spleen, "markers_spleen.csv", row.names = T)

Differentially expressed genes cluster 0 and cluster 3

spl_cl1_vs_cl2 <- FindMarkers(seu_spleen_filt_lymph, ident.1 = 1, ident.2 = 2, only.pos = F)


data.table(spl_cl1_vs_cl2) %>% 
  kbl() %>% 
  kable_styling(font_size = 11,
                bootstrap_options = c("striped", "hover", 
                                      "condensed", "responsive")) %>% 
  scroll_box(width = "100%", height = "400px")

#write.csv(spl_cl1_vs_cl2, "spl_cl1_vs_cl2.csv", row.names = T)

FeaturePlot(seu_spleen_filt_lymph, features = c("Cd3d","Cd4","Cd8a","Sell","Ccl5",
                                       "Foxp3","Il2ra","Pdcd1", "Cxcr5","Klrd1","Klrk1",
                                       "Isg15","Srm","Klrg1", "Mki67"), ncol = 5)

##ggsave("fig/markers1_seu_spleen_filt_lymph.png", width = 20, height = 11)

FeaturePlot(seu_spleen_filt_lymph, features = c("Ifng","Tbx21","Eomes","Bcl6","Ccr8",
                                       "Ccr10","Cxcr3","Tox", "Lag3","Ikzf2","Id2",
                                       "Id3","Prdm1","Itgae", "Cd69"), ncol = 5)

##ggsave("fig/markers2_seu_spleen_filt_lymph.png", width = 20, height = 11)

FeaturePlot(seu_spleen_filt_lymph, features = c("Gata3","Rorc","Il23r","Cxcr6","Cd27",
                                       "Entpd1","Izumo1r","Tnfrsf9", "Tnfrsf18","Tnfrsf14","Il17a",
                                       "Gzma","Gzmb","Gzmk", "Gzmm"), ncol = 5)

##ggsave("fig/markers3_seu_spleen_filt_lymph.png", width = 20, height = 11)

Differential expression

Differentially expressed genes between GFD and STD.

By DeSeq spleen

sum_exp <- AggregateExpression(seu_spleen_filt_lymph, group.by = "MouseID", slot = "counts")

md <- data.frame(Diet = c("GF", "GF", "GF", "GF", "GF", "STD", "STD", "STD", "STD"))

# Create a DESeq dataset
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = sum_exp$RNA,
                              colData = md,
                              design= ~ Diet) 

rownames(dds) <- rownames(sum_exp$RNA)

# Run the DESeq2 algorithm, which will give us the differentially expressed genes
dds <- DESeq(dds)

## PCA
pca = vst(dds, blind = FALSE)
pca_plot <- plotPCA(pca, intgroup = c("Diet"))
pca_plot

## Results data frame
res <- results(dds)
results <- as.data.frame(dplyr::mutate(as.data.frame(res), sig=ifelse(res$padj<0.05, "FDR<0.05", "Not Sig")), row.names=rownames(res))

Genes downregulated in STD

results %>% filter(log2FoldChange < 0) %>%  arrange(pvalue) %>% 
  kbl() %>% 
  kable_styling(font_size = 11,
                bootstrap_options = c("striped", "hover", 
                                      "condensed", "responsive")) %>% 
  scroll_box(width = "100%", height = "400px")